function [] = computeCorrelation_TheoryFix( dirName ,fileNum ,nORr , dataMode)
% compute the correlation function with radius 
% nORr is a parameter can be determined by "plotEntropyFix"
% "dataMode" only be used in function "readText"
%---------------------------Preparation-------------------------------%
% this part repeats the "Preparation" of plotEntropyFix basically
global computeCorrC % has not finished
global topological

[position,velocity,num]=readText(dirName, fileNum,dataMode);

S = sum(velocity) / num; 
SL = norm(S); % SL equals the model of S
unitVec = S/SL; % row vector
sL = velocity * unitVec'; % column vector
pai = velocity - sL * unitVec ;

distance = zeros( num,num );
expCorr = zeros(num,num);
for i= 1 : num
    for j = 1: num
        distance(i,j) = norm(position(i,:)-position(j,:));
        expCorr(i,j) = sum(velocity(i,:).*velocity(j,:));
    end
end % acquire distance and correlation between i and j

[logicI,logicB] = distinguishBoundary(position,distance);
% logicB = zeros(1,num) == ones(1,num);
% logicB(1:100) = true;
% logicI = ~logicB;
numI = sum(logicI);
distance_II = distance(logicI,logicI);

%-----------------------Compute Correlation--------------------%
maxR  =  floor(max(max(distance_II)))+1;
deltaR = maxR / 20;
maxN = floor(maxR/deltaR);

[logic, ~ ,corrSum1] = computeIntCorrelation( num ,nORr ,distance ,expCorr );
N_OO = 0.5 * (logic + logic'); % when i = j ,logic(i,j) = False
N_IB = N_OO(logicI,logicB);
N_II = N_OO(logicI,logicI);
AdJ_II = diag(sum(N_II,2) + N_IB * sL(logicB)) - N_II;
%h_I = N_IB * velocity(logicB,:); % column vector
hL_I = N_IB * sL(logicB);
hP_I = N_IB * pai(logicB,:);
J =(numI-1)*( 0.5 *sum(sum( AdJ_II^-1 .*( hP_I * hP_I'))) - 0.5 * (norm(sum(pai(logicB,:))+sum( AdJ_II^-1 * hP_I)))^2 / (sum(sum(AdJ_II^-1))) +...
    sum(hL_I) + 0.5 * sum(sum(N_OO(logicB,logicB).*(velocity(logicB,:) * (velocity(logicB,:))'))) + 0.5 * sum(sum(N_II))  - 0.5 * sum(corrSum1 .* sum(logic)) )^(-1);

A_II = AdJ_II * J;
hP_I = hP_I * J;

AA = A_II^-1 * hP_I; % controversial expression has been certained
pai_Theo = AA - sum(A_II^-1,2) * ( sum(pai(logicB,:)) + sum(A_II^-1 * hP_I) ) / sum(sum(A_II^-1)) ;

CorrP_Theo = pai_Theo * pai_Theo' + 2 * ( A_II^-1 - sum(A_II^-1,2) * sum(A_II^-1) / sum(sum(A_II^-1)) ) ;
CorrL_Theo = (1-SL-0.5*diag(CorrP_Theo)) * (1-SL-0.5*diag(CorrP_Theo))';
CorrC_Theo = 2 * ( A_II^-1 - sum(A_II^-1,2) * sum(A_II^-1) / sum(sum(A_II^-1)) ); % appears in theory only

CorrP_TheoR =zeros(1,maxN);
CorrL_TheoR =zeros(1,maxN);
CorrC_TheoR =zeros(1,maxN);
count = zeros(1,maxN);

for i = 1 : numI
    for j = 1 : numI
        temp = floor((distance_II(i,j)/deltaR))+1;
        if temp < maxN
        CorrP_TheoR(temp) = CorrP_TheoR(temp) + CorrP_Theo(i,j);
        CorrL_TheoR(temp) = CorrL_TheoR(temp) + CorrL_Theo(i,j);
        CorrC_TheoR(temp) = CorrC_TheoR(temp) + CorrC_Theo(i,j);
        count(floor((distance_II(i,j)/deltaR))+1) = count(floor((distance_II(i,j)/deltaR))+1)+1;
        end
    end
end % transfer (i,j) to R 's relation

meanCorrP_TheoR = CorrP_TheoR./count;
meanCorrL_TheoR = CorrL_TheoR./count;
meanCorrC_TheoR = CorrC_TheoR./count;
% meanCorrP_TheoR(1) =[]; 
% meanCorrL_TheoR(1) =[]; 
%-------------------------------Plot--------------------------------%
numLabel = 1:10:numI;
rLabel = deltaR*(1:maxN);

figure(2)
subplot(1,2,1)
plot(rLabel,meanCorrL_TheoR,'k') % 'kx'
if topological == 1
    text(maxR/2,max(meanCorrL_TheoR),num2str(nORr))
elseif topological == 0 
    text(maxR/2,max(meanCorrL_TheoR),num2str(nORr))
end
subplot(1,2,2)
plot(rLabel,meanCorrP_TheoR,'k')
 
if computeCorrC == 1
    figure(3)
    subplot(1,2,1)
    plot(numLabel,pai_Theo(numLabel),'kx')
    subplot(1,2,2)
    plot(rLabel,meanCorrC_TheoR,'k')
end
% subplot(1,3,3)
% plot(deltaR*(1:maxN),meanCorr_TheoR,'k')

end



